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We present deterministic algorithms for the simultaneous control of an arbitrary number of quan- 
tum observables. Unlike optimal control approaches based on cost function optimization, quantum 
multiobservable tracking control (MOTC) is capable of tracking predetermined homotopic trajec- 
tories to target expectation values in the space of multiobservables. The convergence of these 
algorithms is facilitated by the favorable critical topology of quantum control landscapes. Fun- 
damental properties of quantum multiobservable control landscapes that underlie the efficiency of 
MOTC, including the multiobservable controllability Gramian, are introduced. The effects of mul- 
tiple control objectives on the structure and complexity of optimal fields are examined. With minor 
modifications, the techniques described herein can be applied to general quantum multiobjective 
control problems. 
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I. INTRODUCTION 

The optimal control of quantum dynamics (QC) is re- 
ceiving increasing attention due to widespread success 
in laboratory experiments and numerical simulations 
across a broad scope of systems. With these promising 
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results it becomes imperative to understand the reasons 
for success and to develop more efficient algorithms that 
can increase objective yields, especially for more com- 
plex objectives. 

Quantum optimal control problems studied to date 
fall into two major classes: 1) control of the expecta- 
tion values of single quantum observables or pure states 
P, [H, H, 0h 2) control of quantum dynamical transfor- 
mations [a, The former has been implemented in 
both simulations and experiments, whereas the latter 
has been approached predominantly through numerical 
studies, due to the expense of quantum process tomog- 
raphy. An important third class of QC problems, which 
lies between the latter two problem types, is the control 
of arbitrary numbers of quantum observables. To date, 
few effective techniques - either experimental or numer- 
ical - have been reported for multiobservable quantum 
control 

The typical experimental approach to controlling sin- 
gle quantum observables is to maximize (or minimize) 
a control objective functional, such as the expectation 
value of an observable operator, using stochastic search 
algorithms. A common technique is to randomly sample 
single observable expectation values at various points 
over the control landscape and use genetic algorithms 
(GA) to update the control field [§]. However, the util- 
ity of such stochastic search algorithms needs careful 
consideration for more complex QC problems like mul- 
tiobservable control. A recent work Q examined the 
application of multiobjective evolutionary algorithms 
(MOEA), which do not rely on a control objective func- 
tional, to a two-observable quantum optimal control 
problem. While the MOEAs displayed promising im- 
provements over GAs, they scale more poorly than ob- 
jective function-based algorithms Moreover, both 
GAs and MOEAs are limited in their efficiency by the 
fact that they do not make use of prior knowledge per- 
taining to the structure of the search landscape. As 
such, there remains a need for the development of both 
numerical and experimental control strategies tailored 
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for the problem of multiobservable quantum control. 

Recently, significant strides have been made towards 
establishing a foundation for the systematic develop- 
ment of efficient optimal control algorithms for more 
complex QC problems based on the observation that 
the landscape traversed by search algorithms in the op- 
timization of quantum controls is not arbitrarily compli- 
cated, but rather possesses an analytical structure origi- 
nating in the geometry of quantum mechanics [llj . Prior 
work established important features of these landscapes, 
in particular their critical topologies [!, [HJ ■ The crit- 
ical points of a control landscape correspond to locally 
optimal solutions to the control problem. The land- 
scapes corresponding to problems f ) and 2) were shown 
to have a set of critical submanifolds of measure zero 
in the search domain, and to be essentially devoid of 
local traps, i.e., the vast majority of local suboptima 
are saddles, facilitating the convergence of local search 
algorithms. 

The underlying monotonicity of quantum control 
landscapes has increased interest in deterministic quan- 
tum control o ptim ization algorithms. A recent experi- 
mental study [1 31 ] demonstrated at least a two-fold im- 
provement in optimization efficiency through the use of 
local gradient algorithms rather than a GA. Such experi- 
mental gradient-based search algorithms, enabled by the 
favorable critical topology of QC landscapes, may in fact 
be essential for high-precision quantum multiobservable 
control. Still, any optimal control theory (OCT) strat- 
egy based on optimizing a control cost functional may 
not be ideal for multiobservable control problems, where 
the most effective combination of observable expectation 
values is not always clear at the outset. A more power- 
ful deterministic algorithm would offer the ability to ex- 
plore arbitrary trajectories across the control landscape 
to identify desirable solutions. 

In order to address this need, we develop in this paper 
the general theory of quantum multiobservable tracking 
control (MOTC), a control strategy that seeks to drive a 
quantum system along homotopic paths to desired mul- 
tiobservable expectation values. This algorithm is mo- 
tivated by the so-called continuation methodology for 
multiobjective optimization [14f , which is a vector opti- 
mization technique based on the principles of differen- 
tial geometry that serves as an alternative to multiob- 
jective stochastic optimization. The paper is organized 
as follows. In Section II, we derive analytical results 
characterizing the gradient flows of quantum multiob- 
servable control cost functionals, examining the factors 
that govern the convergence of these flows, and showing 
that they follow paths (in the unitary group) that are 
highly Hamiltonian-dependent. In Section III, we review 
the theory of unitary propagator tracking control, high- 
lighting its algorithmic advantages compared to scalar 
objective optimization, as well as its stringent require- 
ments for the regularity of control fields. In Sections IV 
and V, we develop the related theory of multiobservable 
tracking control and introduce the MOTC Gramian ma- 



trix, which characterizes the ability of the search algo- 
rithm to move in arbitrary directions in multiobservable 
space from a given point on the control landscape. In 
Section VI, we describe the methods employed in the nu- 
merical implementation of multiobservable tracking and 
then, in Section VII, we present numerical illustrations 
of MOTC in the case of several model systems, examin- 
ing the effect of multiple observable constraints on the 
structure and complexity of the optimal controls. Fi- 
nally, in Section VIII, we draw general conclusions and 
discuss future directions. 



II. QUANTUM MULTIOBSERVABLE 
CONTROL GRADIENT FLOWS 

A generic quantum optimal control problem can be 
written fill ] : 

max *(U(T)) (1) 

e(t) 

where U(T) is an implicit functional of e(t) via the 
Schrodinger equation for the unitary propagator 

^± = -i H (s(t))U(t), U(0) = I N . 

Here H denotes the total Hamiltonian and e(t) is the 
time-dependent control field. Solutions to the optimal 
control problem correspond to ^j^j = for all t € [0, T]. 
The objective function <E> can take various forms. The 
most common form of $ is the expectation value of an 
observable of the system: 

*(U(T)) = Tr(C/(T)p(Q)C/t(T)e) ) (2) 

where p(0) is the initial density matrix of the system 
and O is the Hermitian observable operator whose ex- 
pectation value is to be maximized [151 ]. 

A natural objective function for multiobservable con- 
trol is a positively weighted convex sum of individual 
observable objectives, i.e., 

m 

®m(U) =J2 a k$k(U), a k >0, (3) 
fe=i 

where <f> k (U) = Tr(Up(0)WG k ), fe = l,2,---,m. The 
goal of a multiobjective optimization problem may be to 
maximize the expectation values of all observables, i.e., 

$(U(T)) = {*i(tf(T)), • • • , $ M (U(T))}. 

Alternatively, the goal may be to target observable ex- 
pectation values Xkt k = 1, • • • , m, in which case objec- 
tive function §S§ can be replaced by 

^M(U) = J2 a ^k(U)- X k\ 2 , a k >0. (4) 

k=l 



3 



In this section, we examine the factors that affect the 
efficiency of algorithms that optimize scalar objective 
functions of the form or © , with a specific focus on 
gradient algorithms. We begin by writing expressions 
for the gradients of these functionals. An infinitesimal 
functional change in the Hamiltonian 5H(t) produces an 
infinitesimal change in the dynamical propagator U (T) 
as follows: 



U{T)U\t)5H{t)U{t)dt. (5) 



mT) = -iL 

The corresponding change in <E> is then given by 



5$ = -- 

h 



Tr 



([Q(T),U\t)5H(t)U(t)}p(0) 



dt. 



where 9(T) = W{T)OU{T). In the special case of the 
electric dipole formulation, the Hamiltonian becomes 



H(t) = H - n ■ e(t) 



(6) 



where Hq is the internal Hamiltonian of the system and 
[i is its electric dipole operator. In this situation, the 
gradient of $ is Q: 



5)=-^ [«< 7 ->-"<"l"«" 



(7) 



where fi(t) = W (t) \xU (t) , and the gradient for $m is 

r , rn 
8®m I 



Se{t) 



Y,a k Tr([e k (T),[x(t)]p(0)). (8) 



The flow trajectories followed by these gradients are 
the solutions to the differential equation 



de s (t) 
ds 



<5$ 



(M) 



Se s (t) 



(9) 



where s > is a continuous variable parametrizing the 
algorithmic time evolution of the search trajectory, and 
7 is an arbitrary positive constant that we will set to I. 
Prior work [llj ] demonstrated that under the assump- 
tion of linear independence of the elements of the time- 
evolved dipole operator fi(t), the landscape for objective 
function ^ contains no local maxima, thus ensuring 
that the gradient flow © cannot be trapped. Further- 
more, it was shown that the critical set of this objec- 
tive function consists of submanifolds of Lebesgue mea- 
sure zero in U(N) — indicating that the likelihood of 
encountering a suboptimal critical point is essentially 
null. Since equation ([5]) is the gradient of the expecta- 
tion value of a single observable 0m = 5TJfe=i a k®k, it 
follows that its flow will also converge to the global op- 
timum of the objective function and share the above fa- 
vorable landscape features. These features indicate that 
gradient-based algorithms may also be effective for mul- 
tiobservable control optimization. 



The development of such deterministic search algo- 
rithms for quantum multiobservable control is espe- 
cially important because of the aforementioned difficul- 
ties in sampling multiobjective landscapes with stochas- 
tic techniques. However, there are two characteris- 
tics of the simple gradient flows ([9]) that could be 
improved to render them more efficient in searching 
multiobservable control landscapes. First, the con- 
vergence rate of gradient flow control optimization is 
highly Hamiltonian-dependent. To explicitly isolate the 
Hamiltonian-dependent contribution to the search dy- 
namics, consider first the gradient flow of <&m on the 
unitary group U(N), which is given by [161 ]: 

^ = V$ M (V S ) = jTa k [Q k ,V s p(0)Vj]V s (10) 
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= [G M ,V 3 p(0)v2]V s , 



(11) 



with V £ U(N). V<&m(0 denotes the gradient of the ob- 
jective function on U(N), where the Riemannian metric 
is defined by the inner product 

(X,Y)=Tr(X^Y), 

for any X and Y in the tangent space TyU = {VB \ 
£?t = —B} ofU(N) at V. 

We are interested in the relationship between the 
paths followed by the gradient flow on e(t) and that 
(fl"0|) on U{N). The gradient function on e s (t) is related 
to the gradient on U(N) through 

*{™r>>^g>} (12) 



8e s {t) 



A d$ M s(u s (T)); q 

^ d(U s (T)) pq Se s (t) 



(13) 



Now suppose that we have the gradient flow of e s (t) that 
follows (JHJ) and let U S (T) be the projected trajectory on 
the unitary group IA(N) of system propagators at time 
T, driven by e s {t). The algorithmic time derivative of 
U S (T) is then 



ds J Se s (t) ds 

which, combined with ([9]) and (11) , gives 

d(U s (T))ij 



dt (14) 



ds 



T SJJJJJ)) M _ * d$ M 5(U S (T)); 

(15) 



To obtain expressions for the 



note that in the 



electric dipole formulation, equation ([5]) becomes 



6U(T) = -- U(T)U\t)pU(t)8e(t)dt. (16) 
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Taking the derivative of this expression with respect to 
a functional change in the control field, we get 



6U(T) 
5e(t) 



-U(T)U\t)nU(t) 



-f/(T>(i). 



It is convenient to write equation (|15[) in vector form, 
replacing the N x N matrix U S (T) with the N 2 dimen- 
sional vector u.: 



ds 



dt 



l 5e s {t) 5e s (t) 
where the elements of the matrix F are 



V$ju(u s ) :=F[e s (t)]V$ M (u 
(17) 



v, pi 



{U s {T)» s {t)).U s {t)Ut{T)) dt. 



Ipq 



This relation implies that the variation of the propaga- 
tor in U{N) caused by following the gradient flow in the 
space of control fields is Hamiltonian-dependent, with 
the influence of the Hamiltonian completely contained 
in the A 2 -dimensional positive-semidefinite, symmetric 
matrix F[e s (£)]. We will make further use of this decom- 
position in the next section. 

A second drawback to using gradient flows to search 
quantum multiobservable control landscapes is that 
their convergence rate depends on the properties of the 
observables 0& and the initial density matrix p(0). This 
effect is purely kinematic and does not depend on the 
Hamiltonian. In the Appendix we explicitly integrate 
the kinematic gradient flow (11) in the special case that 
p(0) is a pure state, and show that it follows a con- 
voluted path in U(N). In general, the kinematic flow 
evolves on the interior of a polytope whose dimension 
(and the mean path length of the flow trajectory) rises 
with rank and eigenvalue nondegeneracy in p(0) (Ap- 
pendix A). Moreover, it can be shown (Appendix B) 
that each term in the multiobservable gradient can 
be expanded on a "natural" set of basis functions con- 
sisting of linear combinations of matrix elements of the 
time-evolved dipole operator fi(t); the dimension of this 
basis is: 



D = iV 2 -(/V-n) 2 -^', 



n(2iV- 



i=l 



nl (18) 



where n is the rank of p(0) and rii denotes the degeneracy 
of i-th distinct eigenvalue of p(0) (out of the total r such 
distinct eigenvalues). D also increases with the rank and 
eigenvalue nondegeneracy of p(0). 

Gradient flow control optimization is thus 
Hamiltonian-dependent and decreases in efficiency 
for cases where p(0) and the observables Ok are non- 
degenerate and of high rank. Despite these drawbacks 
to gradient flow sampling of multiobservable control 
landscapes, more sophisticated gradient-based algo- 
rithms may offer a significant advantage over stochastic 
search, due to the favorable critical topology of these 
landscapes. We consider these algorithms below. 



III. UNITARY MATRIX FLOW TRACKING 

The symmetric, positive semidefinite matrix F[e s (i)] 
in equation (fT7|) above indicates that the convergence 
time for local gradient-based OCT algorithms may vary 
greatly as a function of the Hamiltonian of the sys- 
tem. Given the decomposition of the gradient into 
Hamiltonian-dependent and Hamiltonian-independcnt 
parts, the natural question arises as to whether the 
Hamiltonian-dependent part can be suppressed to pro- 
duce an algorithm whose convergence time will be dic- 
tated by that of the unitary gradient flow (or a suitable 
kinematic analog), irrespective of the system Hamilto- 
nian. 

In order for the projected flow from e s (t) onto U S {T) 
to match the path followed by the gradient flow (fit))) , the 
quantity that corresponds to movement in each 

step must satisfy a matrix integral equation: 

In the dipole formulation, this relation becomes 

fi s (t)^^-dt = iWl{T)V^ M {U s {T)), (20) 

where, as in the previous section, p, s (t) = Ut(i)uU s (t). 
This is a Fredholm integral equation of the first kind [17| 
for given e s (t) at s and all t G [0,T]. When $ M 

takes the form in equation ©, we have 

p s (t)-^dt = thY,a k [p(p),Ul(T)Q k U.(T)] . 
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On the basis of eigenstates, the matrix integral equation 
([20]) is written 

(i\Mt)\3)^-dt = th(i\(U s (T))^ u (U s (T))\j). 
as 

(21) 

The relation (l2Tj) is underspecified with respect to e s (t), 
indicating that the integral equation possesses a family 
of solutions. To solve it, we must convert it to an ex- 
plicit initial value problem for e s (t) given Eo(t). We first 
expand the partial derivative as 



de s (t) 
ds 



N 

E 



4 j <*M*)li> + /.(*) 



(22) 



on the basis of functions (i\p, s (t)\j), where the "free func- 
tion" f s {t) contains the additional linearly independent 
degrees of freedom. Inserting this expansion into equa- 
tion (j2~l|) produces 



',9=1 



J2 < 9 / mp s (t)\j){p\Mt)\q)dt 



(i\A s \j)- / f s (t)(i\p s (t)\j)dt, 
Jo 
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where A s = ihU}(T)V<f> M (U s (T)). For more general 
target tracks Q s in U(N), we have A s = ^§f- If we 
denote the Gramian matrix G s as 

r-T 



(Gs) 



(*\fta(t)\j)(p \fi a (t)\q)dt, (23) 



II 



we can solve for the coefficients x 



(% a , ■ • • , x 



in..,, . ni . . . 



x? N ) T as 



f=G;1 ' £ flsit ' )fs{t ' )dt ') ' 

provided that G s is invertible. We then obtain the fol- 
lowing initial value problem for s s (t), the algorithmic 
evolution of the control field along the track: 



de s (t) 



/-(*)- 



+ v( ^ - J% s (t')f s (t')dt^j G-'v^m (24) 

where the operator v(-) vectorizes its matrix argument. 
Each "free function" f s corresponds to a unique algo- 
rithmic step in £■(•); modulating this function allows for 
systematic exploration of the set of functions e s (t) that 
are compatible with the gradient step on U(N) [18l ]. 

Solving this set of iV 4 scalar differential equations re- 
quires that the N 2 x iV 2 matrix Gs defined by (f2"3"| is 
invertible, which is equivalent to the claim that the map 
from between control fields and unitary propagators is 
locally surjective in a sufficiently small neighborhood of 
U s ■ Control fields at which G s is singular correspond to 
so-called singular extremal solutions to the control prob- 
lem, in order to contrast them from regular extremals 
[]~9l | . Even if G s is invertible, it is possible that it is 
nearly singular, resulting in large numerical errors dur- 
ing the solution to the differential equation. It is con- 
venient to assess the nearness to singularity of G s by 
means of its condition number, namely the ratio of its 
largest singular value to its smallest singular value. 

Although the gradient function ^f t \ is always lo- 
cally the direction of fastest increase in the objective 
function at s s (t), the path e s {t) (parameterized by s) 
derived from following this gradient has no universal 
(Hamiltonian-independent) global geometry, since <&m 
is not explicitly a function of e s (t). It is known [151 ] 
that this path will not encounter any traps during the 
search, but beyond this, the geometry can be expected to 
be rugged. Unlike the dynamical gradient flow ([9|), the 
algorithmic flow (|2"4")l follows the gradient flow on U(N). 
This flow respects the geometric formulation of the opti- 
mal control objective function in terms of U 8 (T) rather 
directly in terms of e s (t) . The functions fi s (t) contain 
all relevant information about the quantum dynamics, 
whereas the functions ^^jj- contain complete information 
about the geometry of the kinematic search space. The 
N 2 functions (i\/i s (t)\j) , 1 < i,j < N, are calculated 



during the evaluation of g^u\ ! hence, the computational 
overhead incurred by following this flow corresponds to 
that needed to compute the N 4 elements of G s and in- 
vert the matrix, at each algorithmic time step. 

A multitude of other flows could be substituted 

as 



in the RHS of equation (|2Tj) . Since we are interested in 
global optimality, we should choose a flow that follows 
the shortest possible path from the initial condition to 
a unitary matrix that maximizes the observable expec- 
tation value. It can be shown [l5j that a continuous 
manifold of unitary matrices W maximizes These 
unitary matrices can be determined numerically by stan- 
dard optimization algorithms on the domain of unitary 
propagators [20(. The shortest length path in U(N) be- 
tween the initial guess Uo and an optimal W is then the 
geodesic path, which can be parameterized as 



Q s = U exp(iAs) 



(25) 



with A — —i\og(W' ! Uo), where log denotes the com- 
plex matrix logarithm with eigenvalues chosen to lie on 
the principal branch — n < 8 < n. Thus, if we set 
A s = A, the tracking algorithm will attempt to follow 
the geodesic path 1 . 

Due to the nonlinearity of the integral equation (fTTif . 
errors in tracking will inevitably occur, increasing the 
length of the search trajectory beyond that of the mini- 
mal geodesic path. These tracking errors may be amelio- 
rated through the introduction of stabilization terms in 
equation (|24| (see Section [VJ or higher order functional 
derivatives in equation (|20[) . 



IV. MULTIOBSERVABLE HOMOTOPY 
TRACKING CONTROL 

As a methodology for multiobservable control, unitary 
matrix tracking has an advantage over gradient search in 
that it can directly follow an optimal path in the space of 
unitary propagators, assuming the map e(t) <— > U (T) is 
surjective and the first-order formulation of the tracking 
equation (|24]) is sufficiently accurate. However, it can- 
not be implemented experimentally without expensive 
tomography measurements, and carries a computational 
overhead that scales exponentially with system size. 

We now consider a related, experimentally imple- 
mentable tracking control algorithm — multiobservable 
homotopy tracking control (MOTC) — that seeks to 
drive the expectation values of m observable operators 
0i, ...,0 m along a predetermined path w s of $ to de- 
sired target values. These trajectories may correspond 
to expectation value paths corresponding to the kine- 
matic gradient flow (pTJj) , the geodesic (|25|) , or any other 



1 In the case that the control system evolves on a subgroup of 
U(N), e.g. SU(N), the geodesic on that subgroup can be tracked 
instead. 
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path through multiobservable space. In particular, un- 
like the gradient flows ([3]), MOTC may be used to suc- 
cessively drive the expectation values of individual ob- 
servables to their maxima, while constraining the others 
to fixed values. 

The observables measured at each step can be as- 
sumed to be linearly independent without loss of gen- 
erality. Denote the m scalar functions of algorithmic 
time (expectation value paths for each observable) by 
$s, . . . , $"\ Then, the vector Fredholm integral equa- 
tion for analogous to equation (fill)) is given by 



ds 



T 6& a de s (t) dt= dwl 



o 5e s (t) ds 



ds 



l<i<m. (26) 



where $ 4 S = Tr(C/ s (T) ( o(0)D r j'(r)e i ). As in sectionED to 
solve for the algorithmic flow de ^ that tracks this path, 
it is necessary to expand it on a basis of independent 
functions. In this case, we make this expansion on the 
basis of the independent observable expectation value 
gradients: 



de.{t) 
ds 



'6e s (t) 



fs(t), 



(27) 



where the free function is similar to that in equation 
(|22p . Inserting the expansion (f2"T)) into the resulting gen- 
eralized differential equation, we have 



8<S>\ 



o 6e„(t) Se s (t) 



dt 



ds 



<M*) 



fs(t)dt. 



Defining the m-dimensional vector a s (t) 
(aj(t),... ,a™(t)) T by 



5&S 
5e s (t) 



-TT(ul(T)^(t)\7m(U s (T)p8) 



= -^Tr(p(0)[C/t(T)e,C/ s (T) )Ms (t)]), (29) 



and the MOTC Gramian matrix T as 



(T,)« = / al(t)ai(t)dt, 



we can then solve for the expansion coefficients x 
(xl--- ,xT) T as 



x = r„ 



ds J 



provided that T s is invertible. Returning to the original 
expansion (|27|) . we obtain the following nonsingular al- 
gebraic differential equation for the algorithmic flow of 
the control field: 



de a {t) 
ds 



/.(*)• 



ds 



1 T 



a s (t')/ S (t')dt' 



(30) 



In the special case where only one observable is 
measured at each algorithmic step, this equation reduces 
to: 



^ - j\ s (t')f s (t')dt' 



8j^ = 
ds Is x 

'(31) 

where P s is the desired track for (0(T)), and 7 S = 
J Q a 2 s (t)dt. We note that unitary matrix tracking can 
also be expressed as a special case of equation (|30l) . if 
the functions <&* are taken to be the matrix elements of 
U.(T), i.e., 



4>U-W+k(u a (T)) = V\U,(T)\k), j,k= 1, 



,N. 



Then, according to equation (|28[> . the 7V 2 -dimensional 
complex vector a s (t) = v(fi s (t)) is the vector form of 
the matrix fx s (t): 



(i-l)JV+fc 



(*) = -rO'M*)l*), 3, k = !,-••, N. 



With this choice, the N 2 x iV 2 Gramian matrix T is iden- 
tical to G s as defined by (|2"3"j) . Indeed, with <E>* and a 1 
suitably defined, the MOTC tracking differential equa- 
tion (I30p provides a unified framework for continuation 
approaches [3] to generic quantum multiobjective con- 
trol problems. 

One disadvantage of multiobservable tracking, com- 
pared to local gradient optimization of objective func- 
tion $>m, is that it encounters singular control fields 
more frequently. In this case, singularities correspond 
to the situation that variation of control fields near 
the control under consideration is not sufficient to pro- 
duce arbitrary variation of $({7(T)) driven by this con- 
trol. Nonetheless, the likelihood of the matrix T s be- 
ing ill-conditioned - even at singular extremal control 
fields s(t), where G is singular - diminishes rapidly with 
N 2 — m (Section EH}. 



Auxiliary penalties on e s (t), such as practical experi- 
mental constraints on the total field fluence, act to de- 
crease the degeneracy in the solutions to the above sys- 
tem of equations for dS g^ ■ It can be shown [2l[ that 
the choice: 



/.(*) 



e s (t)w(t), 



(32) 



for the free function in the tracking differential equa- 
tions, where w(t) > is an arbitrary weight function 
and the rj s term (typically constant) controls numerical 
instabilities, will determine the d£ g^ at each algorith- 
mic time step s that minimizes fluence. 



V. ERROR CORRECTION 

Errors can occur when tracking paths on hi {N) or sub- 
spaces of its homogeneous spaces for two reasons. First, 
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the algorithmic steps on these spaces will be first-order 
approximations to the true increments (Q 



w 



S k+1 ^S k U' 

w Sfc ) due to discretization error; this error will 
increase as a function of the curvature of the flow trajec- 
tory Q s or w s at algorithmic time s&. Second, the track- 
ing integral equations are formulated in terms of only 
the first-order functional derivatives S V 3 1? or ; the 
truncation error will depend on the magnitude of higher 
order functional derivatives. 

In numerical simulations, error-correction methods 
can be applied to account for these deviations from the 
track of interest. Generally, we choose an error cor- 
rection term c s that reflects the difference between the 
current value of the tracking vector and its target value, 
such that the tracking differential equation (|30|) becomes 

-i T 



de s (t) 
ds 



/.(*)- 



dw s 
ds 



a s (t')f s (t')dt' 



rr x a.(t) 



(33) 

For unitary matrix tracking, one can follow the 

Us k (T) 

This correction can be im- 



(minimal-length) geodesic from the real point w Sk 



to the track point Q Sk [IE 
plemented through the choice 



Sk+1 



log [Qt k U Sk (T)] 



in the discretized version of (f33|) . 

For general multiobservable expectation value track- 
ing, the vector space within which $ s resides is not a Lie 
group, but rather a subspace of a homogeneous space, 
and consequently it is not as straightforward to apply 
error correction algorithms that exploit the curved ge- 
ometry of the manifold 2 . In this paper, we take as the 
correction term a simple scalar multiple of the difference 
between the actual value and the target track, i.e., 



c s = /3(w s - $ s 



(3 > 0. 



Note that these methods can in principle also be im- 
plemented in an experimental setting. The use of Runge- 
Kutta integrators (RK, see below) to solve the MOTC 
differential equation (|33| with error correction often de- 
creases tracking errors considerably, especially for prob- 
lems involving integration over long s intervals. RK in- 
tegration with error correction is the method of choice 
for numerical implementations of MOTC, although it is 
more difficult to apply experimentally. 



VI. NUMERICAL IMPLEMENTATION 

For illustrations of multiple observable tracking, a set 
of m < N commuting operators 0i, - • ■ , © m was em- 



ployed. 6i was a randomly chosen nondegenerate di- 
agonal matrix, and 02,-- - , O m were sequential pure 
state projection operators in the canonical basis. The 
TO-dimensional vector $ s was constructed by taking the 
trace of the product of each of these observable oper- 
ators with the time-evolved density matrix. Numerical 
solution of the tracking differential equations (f3TH |3"T|) 
was carried out as follows. The electric field e s (t) was 
stored as a p x q matrix, where p and q are the num- 
ber of discretization steps of the algorithmic time pa- 
rameter s and the dynamical time t, respectively (i.e., 
for each algorithmic step Sk, the field was represented 
as a q- vector for the purpose of computations). Start- 
ing from an initial guess e So (t) for the control field, the 
Schrodinger equation was integrated over the interval 
[0, T ] by propagating over each time step, producing 
the local propagator 

U{t j+1 ,tj) = exp [-iH^T/iq - 1) ]. 

The propagation toolkit [23[ was used for this purpose. 
Local propagators were precalculated via diagonaliza- 
tion of the Hamiltonian matrix (at a cost of iV 3 ), expo- 
nentiation of the diagonal elements, and left/right multi- 
plication of the resulting matrix by the matrix of eigen- 
vectors/transpose of the matrix of eigenvectors. This 
approach is generally faster than computing the matrix 
exponential directly. Alternatively, a Runge-Kutta inte- 
grator [24[ can be employed for time propagation. 
The time propagators 

U{t j ,0) = U{t j ,t j - 1 )---U(t 1 ,0) 

computed in step 1 were then used to calculate the time- 
evolved dipole operators fJ.(tj) = U' ! (tj,0)fiU(tj,0), 
which can be represented as a (/-dimensional vector of 
N x N Hermitian matrices. The r Sfc and vector a Sk 
were then computed by time integration of the dipole 
functions with an appropriate choice of function f s (t) 
described above; in the present work, f s was set ei- 
ther to or the expression in ([5^)1 . For tracking of 
unitary flows - either the kinematic gradient flow (11) 
or the geodesic flow - the next point Q Sk on the 
target unitary track was calculated numerically through 
Qs k = <3s fc _i exp(— tA Sk ds), using a matrix exponential 
routine. 

Next, the control field e Sfc (t) was updated to e Sfc+1 (t). 
This step required inversion of the iV 2 x N 2 matrix G Sk 
or to x to matrix r Sfc , which was carried out using LU 
decomposition. The quantities T" 1 , a Sfc and c Sfc were 

used to compute the (/-dimensional vector — ■ One 
of two approaches was used to update the field: (i) a 
simple linear propagation scheme, i.e., 



de s (t) 



ds 



2 A Ricmannian metric on the space of density matrices may be 
defined |2'J , but its application to multiobservable error correc- 
tion when m < TV 2 — 1 is nontrivial. 



or (ii) a fourth- or fifth-order Runge-Kutta integrator. 
The updated control field e Sk+1 (t) was then again used 
to the propagate the Schrodinger equation. 
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Adaptive control of step size was used to accelerate 
convergence for both the gradient flow and MOTC al- 
gorithms. For the gradient flow, the minimum of the 
objective function <I> or <!>m along the direction of the 
gradient at step s was located by first bracketing the 
minimum in that direction by parabolic extrapolation, 
followed by application of Brent's method for inverse 
parabolic interpolation to pinpoint the minimum. For 
both MOTC and the gradient flow, the Fehlberg embed- 
ded Runge-Kutta integrator (ASRK5) with Cash-Karp 
parameters |24j was used, which embeds a fourth-order 
integrator within a fifth order integrator, and compares 
the differences between the fourth- and fifth-order es- 
timates to assess truncation error and adjust As. To 
compare efficiency of MOTC and gradient control algo- 
rithms, the minimum step size and error tolerance level 
in ASRK5 were set to the largest values that permitted 
stable integration of at least one of the algorithms. 



VII. EXAMPLES 

As an illustrative example of multiobservable con- 
trol, we examine the problem of targeting the vector 
of m observable expectation values associated with the 
unitary propagator W to which the kinematic gradi- 
ent of a single observable converges. We seek to track 
the multiobservable path w s of expectation values that 
corresponds to the geodesic between Uq and W, i.e. 

w k s = Tr^p{0)cxp{-iAs)U^e k U o exp{iAs)Y with A = 

—i\og(W^Uo), 1 < k < m for various numbers of ob- 
servables m. Although this is a simple incarnation of 
multiobservable quantum control, it is well-suited for il- 
lustration of its basic principles, since it deals with a 
universally applicable control objective, and allows for a 
systematic study of the effects of the imposition of addi- 
tional observable objectives on the search dynamics and 
optimal control fields. 

The examples below employ an 11-dimensional Hamil- 
tonian of the form ([5]), with 



diag {0.1,0.2, • ■ • , 1.1}, 

1, i = j; 

0.15, \i-j\ = l; 

0.08, \i-j\=2; 

0, otherwise. 



(34) 
(35) 



We first assess the properties of the Gramian matrices 
G and Gamma for such systems. 



A. MOTC Gramian matrix estimation 

An important consideration in any deterministic mul- 
tiobservable quantum control optimization is the singu- 
larity of control fields applied to the quantum system, 
where the singularity of the mapping e(t) t— > U(T) (as- 
sessed through the Gramian G) must be distinguished 



from that of the mapping e(t) i— > $(f/(T)) (assessed 
through the Gramian T). Nonsingularity of the latter 
corresponds to the ability to move, through an infinitesi- 
mal change of the control field s s (t) — > e s+ d s (t), between 
two infinitesimally close vectors of multiple observable 
expectation value vectors, <f> — > $ + d$. The Gramian 
G depends only on Hq, p, and T, whereas T additionally 
depends on the eigenvalue spectra of p(0) and ©/,.. Ex- 
perimentally or numerically, an ill-conditioned Gramian 
matrix will result in large tracking errors for e s (t). 

The requirements for a control field to be regular for 
the mapping e(t) U (T) are, in general, more stringent 
than those for the mapping e(t) i— > $(t/(T)), since mul- 
tiobservable control requires control of a subset of the 
parameters of U(T). However, if p(0) is rank-deficient, 
the multiobservable Gramian matrix T can become more 
ill-conditioned than G, since certain paths U S (T) cannot 
be accessed. Fig. [1] compares the condition number dis- 
tributions for G and T for various p(0), where the num- 
ber of observables m = 10, for randomly sampled fields 
e(i) of the form 



N 



N 



e(t)=X! Aijsiniuiijt + faj), < t < T, (36) 
i=i j=i+i 

where w, 3 = \Ei — Ej \ denote the transition frequencies 
between energy levels E i: Ej of H , fa* denotes a phase 
sampled uniformly within the range (0, 27r], and Aij de- 
notes a mode amplitude sampled uniformly within the 
range (0, 1]. The final time T was chosen to be suffi- 
ciently large to achieve full controllability over U{N) at 
i = T[ll|. 

The absolute magnitudes of the condition numbers 
for a given p(0) spectrum depends on the Hamiltonian, 
with random, dense Hq, p, matrices generally exhibit- 
ing more well-conditioned Gramian matrix distributions. 
Typically, a Gramian condition number C > 10 8 results 
in large numerical errors upon inversion, and would be 
expected to compromise the accuracy of tracking (due 
to the sparseness of control field increments Se(t) that 
are capable of driving the system to the corresponding 
state). 

For small numbers of observables m, rank-deficiency 
in p(0) does not shift the condition number distribu- 
tion for G toward substantially higher values, compared 
to that of r, since the number of parameters of p(T) 
that must be controlled is small. For example, the pure 
and thermal mixed states in Fig. [1] both have well- 
conditioned Gramians T, which permit accurate mul- 
tiobservable tracking (see below). 

For a quantum ensemble in thermal equilibrium with a 
bath at temperature T e , the eigenvalues are determined 
by the Boltzmann distribution, i.e., 



A. = 



exp(-^/fcT e ) 

E?=ieM-E l /kT e y 



,N. 



and MOTC can be carried out without the additional 
overhead of density matrix estimation. Through its ef- 
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Fig. 1: Gramian matrix condition number distributions for the 11-level system in section [VIII The amplitudes and phases 
of the modes of the control field e(t) were sampled randomly from the uniform distributions (0, 1] and (0, 2n], respectively, 
with the mode frequencies tuned to the transition frequencies of the system. (A) MOTC Y matrix with 10 observables and 
p(0) being full-rank and nondegenerate (thermal); (B) MOTC V matrix with 10 observables and p(0) being a pure state; (C) 
Unitary tracking G matrix. 



feet on the eigenvalues of p(0), the spectrum of Hq (in 
particular, the energy-level spacings) plays an impor- 
tant additional role in determining the Gramian matrix 
condition number distribution for thermally prepared 
states. 

The regularity of control fields generally becomes 
harder to achieve for larger Hilbert space dimension. 
For example, for Hamiltonians of the form ([6]), the mean 
G condition number rose from O(10 18 ) for N = 11 to 
O(10 19 ) for N = 19, and the T condition number for full 
rank p rose from O(10 5 ) to O(10 6 ) for tracking of one 
complete measurement basis (m — 10, m = 18, respec- 
tively) . 

The distributions above provide rough estimates for 
probabilities of encountering singularities during the im- 
plementation of the different forms of MOTC examined 
in the following sections. However, the actual changes 
in r that occur along an optimization trajectory will 
not be identical to those obtained by random sampling 



from these distributions, since the frequencies u> of suc- 
cessive control fields in MOTC are not independent ran- 
dom variables 3 As such, the fields £ s (t) sampled dur- 
ing MOTC tracking will typically be more regular than 
those in the distributions displayed above, although the 
trends will be similar. 



B. Multiobservable tracking 

According to equations ([7]) and (fl~8"|) , the efficiency of 
gradient flow-based observable maximization decreases 



3 The maximum change in T (and hence its condition number C) 
that can occur along an optimization trajectory is bounded due 
to the fact that the norm of the gradient jj^fy is bounded by 
the norm of the dipole operator Q|. 
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Fig. 2: Comparison of MOTC of 2 (dotted), 4 (dashed) and 
10 (solid) observables from a single measurement basis, for 
11-level system (I34|l . p(0) was of rank 7 with nondegener- 
ate eigenvalues. The target was a unitary propagator that 
maximized the expectation value of the first observable. A) 
Total distance (Frobenius norm) in U(N) traversed by the 
dynamical propagator U(T) versus algorithmic step; B) Dis- 
tance between target Q s and actual U„(T) tracks as function 
of algorithmic step; C) Total distance (Euclidean norm) in 
L 2 (R) between the current field e s (-) and the original field 
Eo(0 along the optimization trajectory. 



with higher rank and nondegeneracy in p(0). In these 
cases, the gradient flow typically follows a longer path 
in both U(N) and e(t). On the other hand, the uni- 
tary path can be constrained more effectively by MOTC 
in such cases, with the most stringent control possible 
for full-rank, nondegenerate p(0). Fig. [5] compares the 
pathlcngths in U(N) of MOTC optimization trajecto- 
ries in such a case with 2, 4 and 10 observable expecta- 
tion values tracked along the geodesic to an optimal W 
precomputed by following the kinematic gradient flow 



100 




0.2 0.3 
Frequency 

Fig. 3: Comparison of optimal fields obtained via m — 2 
(dashed) and m = 10 (solid) MOTC search algorithms that 
maximize a single observable O. p(0) was a thermal mixed 
state. A) Optimal control fields for 2 and 10 observable 
tracking; B) Fourier power spectra of the optimal control 
fields. Note the increase in high frequency modes due to the 
imposition of additional observable objectives. eo(t) was of 
the form (l36l). 



= [9i,V.p(0)V2]V a 4 . The Hamiltonian was 
employed for these simulations. Across all cases, the 
Gramian T s was well-conditioned at each step of MOTC 
optimization, as predicted by the T condition number 
distributions above. 

As can be seen from the Figure, the pathlength in 
IA(N) decreases progressively with increasing to; this 
pathlength is in all cases smaller than that of the gra- 
dient flow (data not shown). For small to, there are 
significant stochastic fluctuations in the unitary step- 
size per iteration, which are smoothed out for larger to. 
By contrast, the Euclidean pathlength traversed in the 
space of control fields e(-) (assessed in terms of the Eu- 



4 Note that although the tracked unitary path Q s is signifi- 
cantly shorter than that followed by the dynamical gradient, still 
shorter paths could be identified by minimizing the distance of 
W to Uq while constraining (0) to remain at its maximal value. 
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clidean distance |e sl (-)-e S2 (-)| 2 = J2i l £ si(*i)- £ s 2 (*i)| 2 j 
where Si, S2 label successive points along the optimiza- 
tion trajectory) increases systematically with m, with 
the smallest change in e(-) along the optimization trajec- 
tory occurring for the gradient and the greatest change 
occurring for m = 10. Almost universally, imposing ad- 
ditional observable tracks increases the field pathlength 
and distance between £o(') and e s (-)- Since the MOTC 
field update (|5D|) is based on projections the multiob- 
servable gradient ((8|) , each additional observable tracked 
causes the optimization trajectory to deviate further 
from that followed by the gradient. 

Fig. [3] depicts the optimal control fields identified by 
m = 2 and m = 10 MOTC algorithms. Note that, even 
in cases such as this where the imposition of multiple 
observable expectation value constraints drives the sys- 
tem to a final unitary propagator W that is closer to 
Uq, the optimal fields are invariably more complex, with 
the Fourier spectra displaying higher frequency modes 
for additional observables. The dissimilarity of these 
fields from the initial guess indicates that optimal fields 
for multiobservable quantum control cannot be identi- 
fied by heuristic reasoning based solely on, for example, 
the spectrum of the Hamiltonian. 

One disadvantage of tracking-based optimization of 
observable expectation values, compared to steepest as- 
cent search, is that more precise measurements of the 
gradient are required to remain on the desired track. 
Tracking paths for additional observables, where the 
auxiliary tracks are chosen so as to constrain the uni- 
tary propagator to more a uniform path, can confer ad- 
ditional stability and robustness to observable tracking- 
based optimization. 

Errors in tracking that occur due to breakdown of the 
first-order MOTC approximation may be stabilized by 
the use of additional observables because it is less likely 
that the first-order approximation will break down si- 
multaneously for all observables. Since more auxiliary 
observables can be tracked when p(0) has more inde- 
pendent parameters, this stabilization method is most 
effective when p(0) is full-rank. Thus, in multiobserv- 
able tracking, a greater rank of p(0) can increase the 
the maximal stability of the algorithm, as well as the 
diversity of possible multiobservable expectation value 
targets and the freedom to follow arbitrary paths to- 
ward those targets. Fig. 0] compares the tracking errors 
for the MOTC variants depicted in Fig. [2] 

The relative efficiency of MOTC and gradient flow 
observable optimization can be assessed by employing 
ODE integrators with adaptive stepsize. For this pur- 
pose, the Fehlberg embedded Runge-Kutta (ASRK5) 
method was used. Although it is somewhat more dif- 
ficult to implement RK integrators experimentally, they 
provide a measure of the maximum stepsize that can 
be taken along a given optimization trajectory without 
incurring errors above a specified tolerance. Fig. [5] com- 
pares the number of ASRK5 iterations needed to solve 
the above maximization problem by integrating a) the 
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Fig. 4: Observable tracking errors for the control problems 
displayed in Fig. [2] over the optimization trajectory. Dotted 
= 2 observables; dashed = 4 observables; solid = 10 observ- 
ables. 

gradient flow ODE © and b) the MOTC flow §U§ with 
w s set to the 10-observable expectation value track cor- 
responding to the geodesic between Uq and W. Another 
adaptive step-size algorithm for improving optimization 
efficiency, a line search ( Section IVij) . was also used in 
the case of the gradient (data not shown; whereas ei- 
ther ODE can be integrated via ASRK5, only the gra- 
dient flow can be simply implemented via a line search 
routine) . The substantial decrease in the number of re- 
quired iterations in the case of tracking algorithms, in- 
dicative of enhanced global optimality of the path they 
follow in L 2 (M.), is displayed in the Figure. 



VIII. DISCUSSION 

We have presented a class of deterministic algorithms 
for the optimal control of multiple quantum observ- 
ables, based on tracking multiobservable expectation 
value paths to desired endpoints. These algorithms 
leverage the previously reported simple critical topol- 
ogy of quantum control landscapes [ll[ , which indicates 
that gradient-based algorithms will rarely encounter sin- 
gularities. An additional feature of multiobservable con- 
trol landscapes that governs the ability of local search 
algorithms to follow arbitrary paths in multiobservable 
space, the MOTC Gramian matrix, has been identified 
and its properties investigated numerically for selected 
problems. Error correction methods have been described 
that should facilitate the experimental implementation 
of these algorithms in the presence of noise. Moreover, 
a general MOTC framework has been presented that 
extends beyond the specific problem to multiobservable 
control to encompass more general quantum multiob- 
jective optimization problems. Extensions to problems 
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Fig. 5: Algorithmic efficiency of observable control algo- 
rithms assessed through variation of step size. Top: adaptive 
stepsize Runge-Kutta (ASRK5) MOTC integration (10 ob- 
servables); Bottom: ASRK5 gradient flow integration, both 
for the observable maximization problem in Fig. [2] Expec- 
tation value refers to the first observable. 



involving objectives in multiple quantum systems [8( are 
straightforward . 

The performance of MOTC algorithms has been com- 
pared to that of local gradient flow algorithms based on 
scalar objective functions. Although the e(i)-gradient 
flow is always the locally optimal path, its projected 
path in U{N) is generally much longer than those that 
can be tracked by global MOTC algorithms. Even for 
single observable control problems, the latter often re- 
quire fewer iterations for convergence. 

Features of the optimal control fields have been com- 
pared for various numbers of controlled observables. Ap- 
plication of multiple observable expectation value con- 
straints has been shown to generally increase the com- 
plexity of optimal control fields on the level set of a quan- 
tum control landscape [ll| through the introduction of 
higher frequency field modes beyond those present in the 
transition frequency spectrum of the internal Hamilto- 
nian. 

Two MOTC applications of particular importance are 
A) the simultaneous maximization of the expectation 
values of a set of observables and B) the preparation of 
arbitrary mixed states. (A) corresponds to the problem 
of identifying Pareto optimal solutions to the multiob- 
servable control problem, i.e. control fields e{t) such that 
all other fields have a lower value for at least one of the 
objective functions or else have the same value for 
all objectives. Whereas in generic multiobjective prob- 
lems, the Pareto frontier is very difficult to sample due 



to its irregular structure, the simple landscape topol- 
ogy and geometry of quantum optimal control problems 
enables highly efficient methods for Pareto front explo- 
ration based on MOTC. 

Problem (B) requires the control of a large number of 
quantum observables (up to A^ 2 — 1 , for Hilbert space di- 
mension N). For such tasks, the eigenvalue spectrum of 
the initial quantum state plays an important role in de- 
termining the maximum possible optimization efficiency. 
In these cases, it is usually possible to accelerate MOTC 
by choosing an optimal set of observable operators to be 
measured at each step. In this paper, we have focused on 
developing the formalism of MOTC and comparing the 
algorithmic properties governing the efficiency of MOTC 
with those of gradient flow and unitary matrix track- 
ing algorithms. In a follow-up work, we will examine 
strategies for the efficient experimental implementation 
of MOTC, as well as methods for combining quantum 
state estimation with MOTC for effective Pareto front 
sampling. 



APPENDIX A: INTEGRATION OF KINEMATIC 
GRADIENT FLOWS: p(0) PURE 

In this section, we analyze the trajectories followed 
by the quantum observable maximization gradient flow 
(equation (11) in the main text) in order to shed light 
on the factors affecting its convergence rate. This ex- 
pression is cubic in U; however, through the change of 
variables from U S (T) to p s (T) = U„(T)p(0)U$(T), we 
can reexpress it as a quadratic function: 

= -U s (T) P (0)— _„( )DJ(T) 

= p 2 s (T)0~2p s (T)e Ps (T) + epl(T) 

= [ Ps (n[p s (T),e]] 

where s denotes the algorithmic time variable of the 
gradient flow in U(N). This quadratic expression for 
the gradient flow is in so-called double bracket form 

Here, we provide an explicit formula pertaining to the 
analytical solution to the above gradient flow for what 
is perhaps the most common objective in quantum opti- 
mal control theory and experiments, namely the maxi- 
mization of the expectation value of an arbitrary observ- 
able starting from a pure state \i). Using the notation 
\ip s (T)) — U 3 (T)\i), the double bracket flow can in this 
case be written: 

*P = ^p-\i) = I© - (A(T)\e\A(T))i] \A(T)). 

If we define x(s) = (|ci(s)| 2 , . . . , cjv(s)| 2 ), where 
ci (s), . . . , cat(s) are the coordinates of \ip{s)) in the basis 
that diagonalizes O, it can be verified that the integrated 
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gradient flow can be written: 

° 2se -(l Cl (o)|V 



x(s) = 



cn(0)\ 2 ) 



Ef=i I c fc (0)hc 



2 e 2sX k 
2sA» 



(e^| Cl (0)| 2 ,...,e^«| Cjv (0)n 



E*=iK(0)| 2 e 2 ^ 

where Ai, . . . , Ajv denote the eigenvalues of <d. 

The optimal solution to the search problem corre- 
sponds to the basis vector ej„ where has its maximal 
eigenvalue (for example, in the case of population trans- 
fer to a pure state |/), ej* = |/)). The distance of the 
search trajectory to the global optimum of the objective 
(framed on the homogeneous space) can be expressed: 



Ha) 



■f=\\x(s) 



(e 2se a;(0)) 

The time derivative of this distance function is 



dt 



\x(s)-e j 4' z 



''i Jt (0)Efcift. - Xk)e 2sXk x k (0) 



Solving for the zeroes of this time derivative reveals 
that the distance between the current point on the search 
trajectory and the solution can alternately increase and 
decrease with time. Moreover, the same holds for the 
distance to the suboptimal critical points of the objective 
function. As the rank and nondegeneracies in the eigen- 
value spectrum of increase, the number of these sub- 
optimal attractors of the search trajectory increases |l6[, 
resulting in an increase in the algorithmic path length. 



APPENDIX B: NATURAL BASIS FUNCTIONS 
FOR THE DYNAMICAL GRADIENT 



The expression for the dynamical observable maxi- 
mization gradient, Eq.(7) in the main text, can be ex- 
panded as: 



<5$ 



= --Tr{[6(T), M (tM0)} 



-{i\»(t)\j)(j\e(T)\i) , 



i=i j=i 



where the initial density matrix is given in terms of 
its nonzero eigenvalues pi as p(0) = E"=i-P«N)(*I' Pi — 
• • • > Pn > Oj J27=iPi = 1' wnere n denotes the rank 
of jo(0). The assumption of local surjectivity of the map 
e(t) i— » U(T) implies that the functions (i\n(t)\j) are N 2 
linearly independent functions of time. The functions 



(<|e(T)|j>(7lM*)|i)-(<Kt)|j)0'|e(T)|<> 



therefore constitute natural basis functions for the gra- 
dient on the domain s(t). From this the expression (18) 
in the main text for the dimension of the natural basis 
follows. 
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